class verletEspresso
{
private:
	caja *box;
	///Paso de tiempo Dt
	double Dt;
	///Temperatura
	double temp;
	///Coeficiente Langevin gamma
	double gamma;
	///Para el termostato Langevin
	double langevin_pref1, langevin_pref2, langevin_pref2_buffer ;
	double factor;

public:

    verletEspresso(caja *bx, double xgamma, double xtemp, double dt)
    {
        box=bx;
        gamma=xgamma;
        temp=xtemp;
        Dt=dt;
        thermo_init_langevin();
        factor=0.5*Dt*Dt;
    }

    void integrate_vv();
    void thermo_init_langevin() ;
    void thermo_heat_up();
    void thermo_cool_down();
    void rescale_forces();
    void friction_thermo_langevin();
    void rescale_forces_propagate_vel();
    void propagate_vel_pos();

};

inline void verletEspresso::integrate_vv()
{
    /* Prepare the Integrator */
    double vvv=box->part[0].vel[0];
    //cout<<endl;
    if(box->nacc==0)
    {
        thermo_init_langevin();
        thermo_heat_up();
        friction_thermo_langevin();
        #ifdef CONFINED
        interaccionConfined(box);
        #else
        interaccionSimple(box);
        #endif
        thermo_cool_down();
        rescale_forces();
    }

  /* Integration loop */

    /* Integration Steps: Step 1 and 2 of Velocity Verlet scheme:
       v(t+0.5*dt) = v(t) + 0.5*dt * f(t)
       p(t + dt)   = p(t) + dt * v(t+0.5*dt)
       NOTE 1: Prefactors do not occur in formulas since we use
               rescaled forces and velocities.
       NOTE 2: Depending on the integration method Step 1 and Step 2
               cannot be combined for the translation.
    */

    propagate_vel_pos();
    friction_thermo_langevin();
    #ifdef CONFINED
    interaccionConfined(box);
    #else
    interaccionSimple(box);
    #endif
//    force_calc();
    /* Integration Step: Step 4 of Velocity Verlet scheme:
       v(t+dt) = v(t+0.5*dt) + 0.5*dt * f(t+dt) */
    rescale_forces_propagate_vel();

   // for(int i=0;i<box->npart;i++) cout<<scientific<<box->part[i].vel[0]<<endl;
}

inline void verletEspresso::propagate_vel_pos()
{
    int i, j;

    for(i = 0; i < box->npart; i++)
    {
        for(j=0; j < 3; j++)
        {
            /* Propagate velocities: v(t+0.5*dt) = v(t) + 0.5*dt * f(t) */
            box->part[i].vel[j] += box->part[i].force[j];
            /* Propagate positions (only NVT): p(t + dt)   = p(t) + dt * v(t+0.5*dt) */
            box->part[i].pos[j] += box->part[i].vel[j];
        }

        box->boundary(i);
    }
}



inline void verletEspresso::thermo_cool_down()
{
    langevin_pref2 = langevin_pref2_buffer;
}


inline void verletEspresso::thermo_heat_up()
{
    langevin_pref2_buffer = langevin_pref2;
    langevin_pref2 *= sqrt(3);
}

inline void verletEspresso::rescale_forces_propagate_vel()
{
    int i, j, np, c;
    double scale, mass,parke;
    double svxy, svxz, svyz, svzz;
    svxy=svxz=svyz=svzz=0.0;
    parke=0.0;
    //scale = 0.5 * Dt*Dt;
    scale = factor;

    for(i = 0; i < box->npart; i++)
    {
        mass=box->part[i].mass;
        /* Rescale forces: f_rescaled = 0.5*dt*dt * f_calculated * (1/mass) */
        box->part[i].force[0] *= scale/mass;
        box->part[i].force[1] *= scale/mass;
        box->part[i].force[2] *= scale/mass;
        for(j = 0; j < 3 ; j++)
        {
            /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * f(t+dt) */
            box->part[i].vel[j] += box->part[i].force[j];
            parke+=box->part[i].mass*box->part[i].vel[j]*box->part[i].vel[j];
        }
        #ifdef CALC_STRESS
            svxy+=(box->part[i].mass)*(box->part[i].vel[0])*(box->part[i].vel[1]);
            svxz+=(box->part[i].mass)*(box->part[i].vel[0])*(box->part[i].vel[2]);
            svyz+=(box->part[i].mass)*(box->part[i].vel[1])*(box->part[i].vel[2]);
            svzz+=(box->part[i].mass)*(box->part[i].vel[2])*(box->part[i].vel[2]);
        #endif
    }
    #ifdef CALC_STRESS
		///Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy/box->vol/(2.*Dt*Dt)/box->npart;
		box->svxz=svxz/box->vol/(2.*Dt*Dt)/box->npart;
		box->svyz=svyz/box->vol/(2.*Dt*Dt)/box->npart;
		box->svzz=svzz/box->vol/(2.*Dt*Dt)/box->npart;
		//cout<<scientific<<svxy/box->vol/(2.*Dt*Dt)/box->npart<<"  "<<svzz/box->vol/(2.*Dt*Dt)/box->npart<<endl;
	#endif
    ///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/(2.*Dt*Dt)/box->npart;
	///Calculamos la temperatura de la caja
	box->temp=box->ekin/1.5;
	//cout<<box->ekin<<"  "<<box->temp<<endl;
}

inline void verletEspresso::rescale_forces()
{
    int i;
    double scale, mass;
    scale = factor;
    for(i = 0; i < box->npart; i++)
    {
        mass=box->part[i].mass;
        box->part[i].force[0] *= scale/mass;
        box->part[i].force[1] *= scale/mass;
        box->part[i].force[2] *= scale/mass;

    }
}

/** overwrite the forces of a particle with
    the friction term, i.e. \f$ F_i= -\gamma v_i + \xi_i\f$.
*/
inline void verletEspresso::friction_thermo_langevin()
{
    int i,j;
    double massf;
    for(i=0;i<box->npart;i++){
        massf = sqrt(box->part[i].mass);
        for ( j = 0 ; j < 3 ; j++)
        {
            box->part[i].force[j] = langevin_pref1*box->part[i].vel[j]*box->part[i].mass
                + langevin_pref2*(double(rand())/RAND_MAX-0.5)*massf;

        }


    }
}


inline void verletEspresso::thermo_init_langevin()
{
  langevin_pref1 = -gamma/Dt;
  langevin_pref2 = sqrt(24.0*temp*gamma/Dt);
  //langevin_pref2 = sqrt(24.0*temp*gamma/Dt);
}

